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ABSTRACT 

Lattice Boltzmann method is used to simulate forced convection in oil 
sands (low permeable porous geometries). Fluid flows through a sandstone 
look alike square geometry with left wall of the geometry kept at higher 
temperature compared to other walls. Investigation is carried out to study 
influence of increased temperature on flow properties by obserx’ing the 
variation in velocity and temperature profiles for various permeability and 
porosity values, which were varied to match the geometrical properties of oil 
sands. Boundary conditions and the relaxation parameter are suitably defined 
to achieve convergence for low values of permeability. Simulation was carried 
out at low Reynolds number, which however, can be extended to higher values 
of Reynolds number. 
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1. INTRODUCTION 

Use of elevated temperature to influence flow properties is well known in many fields 
that involve transport phenomena in porous media like petroleum extraction, 
underground water flow etc. Considering its wide applications, flow in porous media 
has been studied extensively over many decades. However, difficulties in practically 
replicating the real geometries and restrictions in numerical analysis have reserved a 
wide scope for improvements in studying flow in porous media. Nevertheless, the 
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recent evolutions in computing technology and numerical tools have assisted in 
accomplishing this goal substantially. This has also helped in the evolution of Lattice 
Boltzmann method (LBM) as an effective numerical technique and has attracted a 
great focus for its ability to simulate flow in porous media. The very advantage of 
LBM is its simplicity in implementing complex boundary conditions. The 
construction of bounce back rules (i.e. boundary conditions) that governs the behavior 
of particles at the boundaries along with the relaxation parameter describes a specific 
problem. Thus, LBM is a much preferred numerical method to study flows in porous 
media compared to the conventional numerical techniques that involve solutions of 
Navier-Stokes equation. 

A number of researchers have contributed in investigating forced convection in 
porous media. Seta et al. (2006) applied LBM to simulate forced convection in porous 
media using Brinkman-Forchheimer equation. The paper also justified that LBM can 
simulate natural convection in porous media in both Darcy and non-Darcy regions in 
an elementary volume scale. Guo and Zhao (2002) presented a LBM model for 
incompressible flows through porous media. Mehrizi et al. (2012) studied forced 
convection heat transfer in a square cavity with inlet and outlet slots and hot obstacle 
with and without porous medium. Grucelski and Pozorski (2013) applied single 
relaxation time variant LBM to the two-dimensional simulations of viscous fluids for 
circular cylinder and a computer-created porous medium. The research in this area is 
extensive. A large amount of this research has focused on flow through porous media 
for theoretical permeability values. Practically, since the permeability of real 
geometries li ke sand, oil sands etc. are very low, i.e. in the range of 10 10 to 10 6 [9], 
certain modification in the numerical methods are necessary. This will include 
applying suitable boundary conditions and the relaxation parameter. The present 
paper investigates forced convection in low permeable geometries like oil sands 
(sandstones), sand etc. which are seen to have permeability and porosity in the range 
of 10 6 to 10 10 and 0.4 to 0.6, respectively. 

2. GOVERNING EQUATIONS AND NUMERICAL METHOD 
2.1. Governing Equations 

The generalized model for incompressible flow in porous media is given by Peng et 
al. (2003) 


V.w = 0 

(1) 

1 VM+ v e V 2 u + F 
at \sj p 

(2) 

where p is the fluid density, u and p are the volume-averaged 
pressure, respectively, f is the total body force given by 

velocity and 

F = -—u-^£=\u\u + sG 
K Jk' 1 

(3) 

where v is the shear viscosity of the fluid and G is the body force. Permeability K 
and the Forcheimer's term f e are related to porosity fas described by (Ergun 1952, 

Seta et al. 2006) 


, 1-75 ,, ^d 2 p 

£ Ue 2 ‘ 150(l-f) 2 

(4) 
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where d p is the solid particle diameter. 


2.2. Lattice Boltzmann Method 

The lattice Boltzmann method can be used to model hydrodynamic or mass transport 
phenomena by describing the particle distribution function fj(x,t) giving the 
probability that a fluid particle with velocity <?, enters the lattice site x at a time 
(Mohamad 2011, Chen and Doolen 1998). The subscript ('represents the number of 
lattice links and / = o corresponds to the particle at rest residing at the center. Present 
investigation is carried out for incompressible fluid flows and a nine -velocity model 
on a two-dimensional lattice (D2Q9). The distribution function governed by 

the lattice Boltzmann equation as (Guo and Zhao, 2002) 


fi(x+ e/dt , t + dt) - fj (. x , t ) = — — ^9 — ft^xd) + 

T (5) 

where dt is the time, r is the relaxation time and is related to the kinematic 
viscosity v by 

u = c s 2 dt^t-^j ( 6 ) 

c s is the sound speed expressed by c s = dx/(S\lt = c/S (cis the particle speed and 
dx i s the lattice spacing). f[ q (x,t) is the equilibrium distribution function for D2Q9 
given by 


fi eq (X,t) = Wip{x,t\ \ + -^(e i .u(x,t))+^-^{e i .u(x,t)) 2 


2 ec, 


u(x,t) 2 

2sc] 


where u(x,t) is the velocity and w t is the weight coefficient. 
The force term F t is given by (Guo and Zhao, 2002) 


(7) 


Fj = Wjd 1 

1 2r 


e i-F uF : (efij - c s I 

9 2 

CZ SC 7 


The fluid velocity u is given by 


c o +y c o + c l| M ;| 

where u, is a temporal velocity given by 


( 8 ) 

(9) 


u t = Yjifi + -^P !G 

i 

The two parameters are given by 




c i 


dt_FE_ 

2 


(10) 


( 11 ) 


2.3. Energy Equation 

The thermal lattice BGK model proposed by Peng et al. (2003) is used to model heat 
transfer 
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gi(x+ ejdt, t + dt) — gj ( x , t ) = 


g i l (x,t)-g i (x,t) 


( 12 ) 


g t (x,t) is the thermal distribution function, r„ the relaxation time and g- q (x,t) is the 
equilibrium distribution function given by (Seta et al. 2006, Zou and He 1997) 


g 0 9 (x,t)=- 


2 pe ii 


(13a) 


8f q (x,t) = f 


3 3 e t .u 9 (e^uf 3 u 2 

2 + 2~^ + 2 r 4 2 r 2 


(i = 1,2, 3,4) 


(13b) 




6 9 {ej.ii) 2 3 u 2 

27 f + 2 ~f 2 ^ 


(i = 5, 6,7,8) 


(13c) 


Reynolds number, Re = , where «„ and L is the characteristic velocity and the 


characteristic length, respectively. 


2.4. Boundary Conditions 

Second order bounce back rule for non-equilibrium distribution function f) is used to 
determine velocity on the four walls. The distribution functions are given by Zou and 
He (1997) 


fr=fr (14) 

where b is the opposite direction of a . 

For energy distribution function g t , second order extrapolation rule is used on the 
left wall and the boundary conditions for all other walls were defined as per method 
introduced by D’Orazio et al (2004). The distribution functions on right wall were 
defined as 


*l=f k- T p) S5=^(l + 3( M+ v)) g& =-^T eq {i+3(u-vj) 

9 9 

611 T ) 

where T p = gQ +g 2 +g 3 +g 4 +g 6 + g 7 and T e< > = ; pl . 

1 + 3 u 

The distribution functions on the top wall were defined as 

g 4 =~I’ e9 (l+3w) g 1 =^-T e ‘>{l+3(u+v)) g H =^-T e H\+3(u-v)) 

9 36 36 

9 9 

where Tp ~ s ^ + 8 \+ 82 + 8 i+ 85 + 8 b 


-67’ 

and T = -. 

1 + 3 u 

The distribution functions on the bottom wall were defined as 

82 ='yT eq (\ +3m) g 5 =-Lr^(l + 3( M + v)) g6 =±T ec <(l+3(u-v)) 

9 36 36 


where 


Tp - £o +Sl + §3+^4+§7 + g9 


(15a) 

(15b) 

(16a) 

(16b) 

(17a) 
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and T eq =-^~ (17b) 

1 + 3 u 

2 . 5 . Problem Description and Numerical Implementation 

Fluid flows through a square geometry containing porous media with the left wall 
kept at a normalized temperature and all other walls adiabatic. Fluid enters the 
geometry with uniform velocity u 0 through the left wall along horizontal direction. 
Permeability and porosity varies between 10 10 to 10 6 and 0.4 to 0.5, respectively, to 
match the geometrical properties of oil sands. Reynolds number Re varies from 0.1 to 
10 and Prandl number Pris taken as 3.8 (of water). The relaxation parameter r is has 
is taken as 1.5. 



(a) 



log(grad p) 
(b) 



log(grad p) 


(c) 

Figure 1 Plot of log (q) versus iog(V/pfor Re=0.1 ands = 0.4 for (a) K = 10 6 (b) K = 1CT 

(c)^ = l(T 10 . 
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3. RESULTS AND DISCUSSION 

The starting point for investigating fluid flow through porous media is the traditional 
modeling approach, typically, Darcy law, relating the pressure gradient Vp to the pore 
fluid flux or Darcy flux q as 

q = -f — W - pg) = -f— W 

UJ ( 18 ) 

where q and p are defined in terms of w and p as q = eu and p = pi?, , respectively. 
Darcy law ensures a linear relationship between iog(g)andiog(V/p. Thus, simulation of 
fluid flow in porous media can be validated by showing that the plot of logtg) versus 
log(Vp) gives a straight line with gradient one. 

While the paper carried out numerical investigation to study the influence of 
elevated temperature and the geometrical parameters on flow properties of fluid, the 
numerical strategy was verified by plotting logL/j and iog(Vp)as shown in Fig 2. For the 
sake of space and invariability in values, only the plots for Re=0.1 are presented. 
Table 1 presents the gradient values of log( 4 )andlog(Vp)and a fine agreement was 
observed between the results and the validating criteria of gradient one. Subsequently, 
simulation was carried out for governing parameters like Reynolds number, porosity, 
permeability etc. on the flow properties of fluid flowing through oil sands. This was 
achieved by varying Re between 0.1 to 10 keeping porosity and the permeability of 
the sandstone fixed, and then varying porosity and permeability to characterize the 
nature of the sandstone. 

Fig 2 presents u-velocity profiles for e = 0.4 at the centerline of the geometry as K 
varies from 10 10 to 10 6 at various Re. The u-velocity tends to become uniform along 
the geometry as flow becomes steady, as shown in Fig 2. The vertical component of 
flow is almost negligible at the order of 10 n . As permeability of the sandstones is 
decreased, the velocity profiles follow the same trend, but with u-velocity decreased 
by same magnitude by which permeability is decreased. This implicates the strong 
influence of permeability flow properties of the fluid in sandstones. The magnitude of 
u-velocity is also influenced by Re, as seen in Fig 2, the velocity increases which is 
due to a stronger flow influenced by heat transfer through the fluid. 

Table 1 Gradient values for iog( ? )andiog(Vp) 


s = 0.4 a = 0.5 




Re 



Re 



0.1 

1 

10 

0.1 

1 

10 

II 

© 

1 

OS 

0.9773 

1.0011 

0.9969 

0.9764 

1.0011 

0.9994 

>3 

II 

o 

1 

oo 

0.9764 

1.0011 

0.9994 

0.9764 

1.0011 

0.9994 

K = 10“ 10 

0.9964 

1.0011 

0.9944 

0.9764 

1.0011 

0.9994 
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0 0.2 0.4 0.6 0.8 1 


y/L 

(b) 



0 0.2 0.4 0.6 0.8 1 

y/L 

(C) 

Figure 2 u-velocity profiles at the centerline for various values of Re at e = 0.4 and (a) 

K = lO^ 6 (b) K = lO- s (c)K = W~ 10 . 
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Figure 3 u-velocity profiles at the centerline for various values of Re at e = 0.5 and (a) 

K = 10~ 6 (b) K = 10 “ 8 (c) A' = 10 _1 °. 
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(a) 




Figure 4 Temperature profiles at various cross section for K = 10 6 and s = 0.4 for (a) 

Re=0.1 (b) Re=l (c) Re=10. 
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(a) 



(b) 



(c) 

Figure 5 Temperature profiles at various cross section for K = 10 6 and e = 0.5 for (a) 

Re=0.1 (b) Re=l (c) Re=10. 
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Since the permeability and Re are small, non-linear linear resistance is neglected, 
which become significant at higher values of K and Re. Fig 2-3 also demonstrates the 
influence off on flow properties. As s is increased from 0.4 to 0.5, the u-velocity 
profiles decrease significantly, naturally, because the void space in the geometry 
increases resulting in a weaker flow through it. This agrees with the rheological 
behavior of the fluid when passing though a medium with bigger void portion. Since 
sand and oil sands are believed to have porosity not more than 0.5 (Bejan and Nield, 
2013), further values off are not considered. The temperature profiles at different 
cross sections, as Re varies from 0.1 to 10, are presented for a: = to -6 and e = 0.4 in Fig 
4. As Re increases from 0.1 to 10, the strength of heat transfer increases, indicating a 
significant increase in heat transfer though the fluid. Fig 5 presents the temperature 
profiles at the same cross section for K = to -6 and s = 0.5 . No significant influence of 
variation in K and e is seen on the temperature profiles, and the plots follow a similar 
pattern as for K = 10 6 , thus, not included. This suggest that the variation in temperature 
profiles is only due to the variation in Re. The temperature profiles follow the same 
trend, both in magnitude and variation. These figures suggest that a better way to 
enhance flow through low permeable geometries can be achieved by elevating 
temperature to influence the flow strength. 


Investigation was carried out to study the influence of elevated temperature and 
geometrical parameters on fluid flow in oil sands. LBM was validated as an 
appropriate numerical tool to simulate forced convection in porous media and applied 
to fluid flows in low permeable geometries. The results significantly depend on 
boundary conditions and the relaxation parameter, which had to be appropriately 
defined to simulate flows in such geometries. The results ascertain that elevated 
temperature plays an important role in influencing flow properties in low permeable 
geometries, which may be of significant application in petroleum extraction and 
groundwater flow etc. The geometric parameters of the porous media has very 
negligible influence on the flow properties compared to external influences, and a 
better way to influence the flow properties would be by influencing the strength of the 
flow by elevated temperature. 
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